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Abstract 



The climate belongs to the class of non-equilibrium forced and dissipative 
systems, for which most results of quasi-equilibrium statistical mechanics, in- 
cluding the fluctuation-dissipation theorem, do not apply. In this paper we show 
for the first time how the Ruelle linear response theory, developed for studying 
rigorously the impact of perturbations on general observables of non-equilibrium 
statistical mechanical systems, can be applied with great success to analyze the 
climatic response to general forcings. The crucial value of the Ruelle theory lies 
in the fact that it allows to compute the response of the system in terms of expec- 
tation values of explicit and computable functions of the phase space averaged 
over the invariant measure of the unperturbed state. We choose as test bed a 
classical version of the Lorenz 96 model, which, in spite of its simplicity, has a 
well-recognized prototypical value as it is a spatially extended one-dimensional 
model and presents the basic ingredients, such as dissipation, advection and the 
presence of an external forcing, of the actual atmosphere. We recapitulate the 
main aspects of the general response theory and propose some new general re- 
sults. We then analyze the frequency dependence of the response of both local 
and global observables to perturbations having localized as well as global spatial 
patterns. We derive analytically several properties of the corresponding sus- 
ceptibilities, such as asymptotic behavior, validity of Kramers-Kronig relations, 
and sum rules, whose main ingredient is the causality principle. We show that 
all the coefficients of the leading asymptotic expansions as well as the integral 
constraints can be written as linear function of parameters that describe the 
unperturbed properties of the system, such as its average energy. Some newly 
obtained empirical closure equations for such parameters allow to define such 
properties as an explicit function of the unperturbed forcing parameter alone 
for a general class of chaotic Lorenz 96 models. We then verify the theoretical 
predictions from the outputs of the simulations up to a high degree of precision. 
The theory is used to explain differences in the response of local and global ob- 
servables, in defining the intensive properties of the system, which do not depend 
on the spatial resolution of the Lorenz 96 model, and in generalizing the con- 
cept of climate sensitivity to all time scales. We also show how to reconstruct 
the linear Green function, which maps perturbations of general time patterns 
into changes in the expectation value of the considered observable for finite as 
well as infinite time. Finally, we propose a simple yet general methodology to 
study general Climate Change problems on virtually any time scale by resorting 
to only few, well selected simulations, and by taking full advantage of ensemble 
methods. The specific case of globally averaged surface temperature response 
to a general pattern of change of the CO2 concentration is discussed. We be- 
lieve that the proposed approach may constitute a mathematically rigorous and 
practically very effective way to approach the problem of climate sensitivity and 
climate change from a radically new perspective. 
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1 Introduction 



A crucial goal in the study of general dynamical and statistical mechanical systems 
is to understand how their statistical properties are altered when we introduce a per- 
turbation related to changes in the external forcing or in the value of some internal 
parameters. The ability to compute the response of the system is of great relevance 
for purely mathematical reasons as well as in many fields of science and technology 

The climate system is an outstanding example of a non-equilibrium, forced and 
dissipative complex system, forced in first instance by spatial differences and temporal 
variability in the net energy flux at the top of the atmosphere. On a macroscopic level, 
as a result of being far from equilibrium, the climate system behaves as an engine, 
driven by the temperature difference between a warm and a cold thermal pool, so that 
the atmospheric and oceanic motions are at the same time the result of the mechanical 
work (then dissipated in a turbulent cascade) produced by the engine, and are processes 
which re-equilibrate the energy balance of the climate system [1, 2, 3, 4]. 

A primary goal of climate science is to understand how the statistical properties - 
mean values, fluctuations, and higher order moments - of the climate system change 
as a result of modulations to some crucial external (e.g. solar irradiance) or internal 
(e.g. atmospheric composition) parameters of the system occurring on various time 
scales. A large class of problems - those involving climate sensitivity, climate variability, 
climate change, climate tipping points - fall into this category. In a system as complex 
and as extended as the climate, where lots of feedbacks are active on a variety of 
spatial and temporal scales, this is in general a very difficult task. The need for 
scientific advance in this direction is outstanding as one considers that even after several 
decades of intense scientific efforts, the accurate evaluation of the climate sensitivity par 
excellence, i.e., the change of the globally averaged surface temperature for doubling of 
CO2 concentration with respect to pre- industrial levels (280 ppm to 560 ppm circa), is 
a tantalizing endeavor, and large uncertainties are still present [5]. 

Such efforts have significant relevance also in the context of the ever-increasing 
attention paid by the scientific community to the quest for reliable metrics to be used for 
the validation of climate models of various degrees of complexity and for the definition 
of strategies aimed at the radical improvement of their performance [6, 7}. The pursuit 
of a quantum leap in climate modelling - which definitely requires new scientific ideas 



4 



rather than just faster supercomputers - is becoming more and more of a key issue in 
the climate community [8]. 

A serious, fundamental difficulty in the analysis of the non equilibrium systems is 
that the fluctuation-dissipation relation [9], cannot be applied [10]. This greatly limits 
the ability of understanding the response of the systems to external perturbations by 
looking at its variability. In the specific case of climate, this can be rephrased by saying 
that climate change signals need not project on the natural modes of climate variability. 
The non-equivalence between free and forced climate fluctuations had been suggested 
by Lorenz [11]. The basic reason for this behavior is that, since the dynamics is forced 
and dissipative, with the asymptotic dynamics taking place in a strange attractor, 
natural fluctuations and forced motions cannot be equivalent. Whereas natural fluc- 
tuations of the system are restricted to the unstable manifold, because, by definition, 
asymptotically there is no dynamics along the stable manifold, external forcings will 
induce motions - of exponentially decaying amplitude - out of the attractor with prob- 
ability one. The fluctuation-dissipation relation can be recovered only if we consider 
perturbations with the somewhat artificial property of being everywhere tangent to the 
unstable manifold or, in a more fundamental way, if we add a stochastic forcing, which 
has the crucial effect of smoothing the invariant measure [12, 13]. Potential links to 
these issues can be found in recent papers proposing new algorithms for three [14] and 
four [15] dimensional variational data assimilation, where it is shown that the quality 
of the procedure improves if the increment of the variables due to the assimilation is 
performed only along the unstable manifold. 

Recently, Ruelle [10, 16] introduced a rigorous mathematical theory allowing for 
computing analytically, ab initio, the response of a large class of non-equilibrium sys- 
tems to general external perturbations featuring arbitrary time modulation. The cru- 
cial result is that the changes in the expectation value of a physical observable can be 
expressed as a perturbative series in increasing powers of the intensity of the external 
perturbation, where each term of the series can be written as the expectation value of 
some well-defined observable over the unperturbed state. In a previous paper [17] we 
showed that the Ruelle theory is, thanks to this property, formally analogous to usual 
Kubo response theory [18], which applies for quasi-equilibrium system. The crucial 
difference lies on the mathematical properties of the invariant measure, which is ab- 
solutely continuous in the quasi-equilibrium case and singular in the non-equilibrium 
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case. 

Ruelle's analysis applies for non-equilibrium steady state systems [19] possessing 
a Sinai-Ruelle-Bowen (SRB) invariant measure, often referred to as Axiom A system 
[20, 21]. This class of systems, even if mathematically non-generic, includes on the 
other hand excellent models for general physical systems, as made clear by the chaotic 
hypothesis [22, 23], which can be interpreted as an extension of the ergodic hypothesis 
to non-Hamiltonian systems [19]. See [24] for an original geophysical perspective. 

The Ruelle response theory, with the support of chaotic hypothesis, has interest- 
ing conceptual implications for climate studies. In fact, the possibility of defining a 
response function basically poses the problem of climate response to forcings and of 
climate change in a well-defined context, and, when considering the procedures aimed 
at improving climate models, justifies rigorously the procedures of tuning and adjust- 
ing the free parameters. Moreover, the response theory allows to compute the climate 
sensitivity, in the special case when static perturbations to the system parameters are 
considered. 

Previously, a response formula was proposed by Cacuci for evaluating the linearized 
change of the solution of a time-independent generic system of nonlinear equations as 
a result of a change in the system's parameters [26, 27]. This can be interpreted as 
a special case of Ruelle's theory, where the unperturbed attractor is constituted by a 
fixed point and a static perturbation to the system evolution equation is considered. 
Cacuci proposed to study this problem using the adjoint operator to the original sys- 
tem, which provided an efficient way to determine the impact of small perturbations. 
Interestingly, early prominent applications of the so-called adjoint method and its ex- 
tension to time-dependent problems, which allowed for evaluating all possible linear 
sensitivities of an evolving model in just one simulation, were been proposed for cli- 
mate related problems. In particular, it was used to evaluate the sensitivities of a 
simple radiative-convective model possessing an attractor constituted by just one fixed 
point [28, 29], and later, in an inherently heuristic way, for studying the response of 
a (chaotic) simplified general circulation model to doubling of the CO2 concentration 
[30]. Whereas the adjoint method did not find much space in further climatic studies, 
mostly due to early discouragement for the computational burden of constructing the 
suitable operators for evaluating the sensitivities, it subsequently reached great success 
in data assimilation problems for geophysical fluid dynamics [31, 32], to the point that 
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a tangent and adjoint model compiler able to automatically generate adjoint model 
code was has been introduced [33]. More recently, a link between advanced adjoint 
techniques and the Ruelle theory has been proposed [34]. 

In the last decade on one side a great effort has been directed at extending the Ruelle 
response theory for more general classes of dynamical systems (see, e.g., [35, 36]), and 
recent studies [17] have shown that, thanks only to the causal nature of the response, it 
is possible to apply all the machinery of the theory of Kramers-Kronig (KK) relations 
[37, 38, 39, 40] for linear and nonlinear processes to study accurately and rigorously the 
susceptibilities describing in the frequency domain the response of a general observable 
to a general perturbation. 

Moreover, the actual applicability of the theory has been successfully tested in a 
number of simple dynamical systems case for the linear [41, 42] and nonlinear [43] 
response. Such numerical investigations have clarified that even in systems which are 
not Axiom A, like the Lorenz 63 system [44], it is possible to successfully use the 
response theory to construct linear [41] and nonlinear susceptibilities [43] which obey 
all of the constraints imposed by the KK theory up to a high degree of precision. 

These investigations definitely motivate further studies aimed at understanding to 
what extent the response theory is an efficient tool for analyzing complex and com- 
plicated systems. In this paper, we take up such a challenge and consider the Lorenz 
96 (L96) system [45, 46, 47], which provides an excellent and celebrated prototypical 
model of a one dimensional atmosphere. The variables of the L96 model can be thought 
as generic meteorological quantities extending around a latitudinal circle and sampled 
at a regular interval. In spite of not being realistic in the usual sense, the L96 model 
presents the basic ingredients, such as dissipation, advection and the presence of an 
external forcing, of the actual atmosphere. For this reason, L96 has quickly become 
the standard model to be used for predictability studies [48, 49, 50], when testing data 
assimilation techniques [14, 15, 51], and new parameterizations [52]. 

Although we are unable to prove that the unperturbed L96 is an Axiom A system, 
in general and for the specific choice of parameters used in our numerical simulations in 
particular, we adopt the chaotic hypothesis and present the first thorough investigation 
of a spatially extended system by using the rigorous statistical mechanical methodolo- 
gies presented in [10, 16, 17, 43]. Moreover, since L96 is a spatially extended system, 
we also explore the applicability of the response theory in all possible combinations 
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of global/local observables and global/local perturbations. We compute rigorously the 
corresponding linear susceptibilities, verify the KK relations and the related sum rules, 
and find an empirical power law, which, as in [53], supports the validity of the chaotic 
hypothesis, allowing to extend the results obtained for our specific choice of model's 
parameters to a rather general class of L96 systems. We also show how to go from the 
frequency back to the time domain, thus deriving from the susceptibility the Green 
function, which acts as time propagator of the considered perturbation for the con- 
sidered observable. The Green function allows to predict, in an ensemble mean sense, 
the change in the observable at any time horizon as a result of a perturbation with 
the same spatial patter as that considered in the calculation of the susceptibility but 
featuring a general time modulation. 

Finally, building upon the results presented here, we propose a simple yet general 
methodology to study general Climate Change problems on virtually any time scale 
by resorting to only few, well selected simulations, and by taking full advantage of 
ensemble methods. The specific case of globally averaged surface temperature response 
to a general pattern of change of the CO2 concentration is discussed. 

Whereas the paper aims at proposing new methods for tackling classical problems 
of climate science, most of the results and of the methodologies proposed are of more 
general interest. In this paper we limit our attention to the linear response. We refer to 
[17, 43] for a theoretical and numerical studies of higher-order effects of perturbations. 

The paper is organized as follows. In Sec. 2 we briefly analyze the general theo- 
retical background of the linear response theory and of the properties of the frequency 
dependent susceptibility and present some new useful results. In Sec. 3 we present the 
main features of the L96 system, introduce the considered perturbations to the forcing, 
derive some basic properties of the response of various observables, and present the 
theoretical predictions. In Sec. 4 we present the results of our numerical investigations 
and describe how they can be generalized to the entire family of L96 models. In Sec. 
5 we provide a relevant example to illustrate how the results presented in this paper 
can be used to devise simple yet rigorous methods to study the climate response at all 
time scales on models of any degree of complexity. In Sec. 6 we discuss the conclusions 
and present perspectives for future work. 
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2 Theoretical Background: Ruelle Theory and Dis- 
persion Relations 



2.1 Definition of the Linear Susceptibility 

We consider an Axiom A dynamical system described by the evolution equation x = 
F(x), so that the invariant probability measure po of the associated flow is of the SRB 
type [10]. Let ($) be the expectation value of the general observable <3> defined as 
J po{dx) We perturb the flow of the system by adding a on the right hand side 

of the evolution equation a vector field X(x)f(t), where X(x) defines the pattern of 
the perturbation, and f(t) is its time modulation. The resulting evolution equation 
results to be x = F(x) + X(x)f(t). Following Ruelle [10], we express the expectation 
value of $(cc) in the perturbed system using a perturbative expansions as: 

oo 
n=l 

Each term of the perturbative series can be expressed as an n-convolution integral of 
the n th order causal Green function with n delayed time perturbation functions [25, 17]. 
Limiting our attention at the linear case we have: 

/+oo 
daiGSVO/ft-'i) (2) 
-oo 

The first order Green function can be expressed as follows: 

= f Po (dx)e(a 1 )AU (a 1 )^(x), (3) 

where A(«) = X(x)V(«) takes into account the effects of the perturbative vector field, 
is the usual Heaviside distribution and n the unperturbed time evolution operator 
so that UqF(x) = F(x(t)) for any function F, with x(t) following the unperturbed flow. 
Note that it is possible to express the Green function as the expectation value of a non- 
trivial but computable observable over the unperturbed SRB measure po- Therefore 
the knowledge of the unperturbed features of the flow is sufficient to define the effects of 
any external perturbation over any observable of our system. In the frequency domain 
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we find that the first term of the perturbative series can be written as: 




fa>iX* MfM x 5(u - wi) = X* H/(w), 



(4) 



where the Dirac delta implies that we are analyzing the impact of perturbations in the 
frequency-domain at the frequency oo. The linear susceptibility is defined as: 



It is important to underline with a thought experiment the computational relevance of 
the last equations and the importance of the susceptibility function. 

Let's suppose we introduce a time dependent perturbation f a (t) to a given pattern 
of forcing X(x), simulate the system and observe the time response of an arbitrary 
observable We now compute the Fourier transform of the observed signal 

and of the forcing modulation. Inverting Eq.(4), we can find the linear susceptibility 

(oo). Let's now consider a different time-modulating function of the forcing fp(t) and 
its corresponding Fourier transform fp(oo). Taking into account Eq. (4), if we multiply 
fp(oo) times the previously computed function x^i^) we directly obtain (^p)^(oo), 
the frequency-dependent response of the observable $ to the forcing X(x) modulated 
by the new function. By applying the inverse Fourier transform we obtain the time- 
dependent response (<f> p)^ (t) without needing any additional simulation. 

Moreover, the knowledge of the susceptibility function allows us to reconstruct the 
G$(t) by inverting Eq. (5). Otherwise, the Green function can be obtained directly 
from observing the response signal by performing a simulation where f(t) = 5(t): in 
this case (see Eq. (2)) we simply have (i) = G^\t). 

2.2 Kramers-Kronig relations and sum rules 

As we see from Eq. (3), for an arbitrary choice of the observable and of the perturbation 
the corresponding linear Green function is causal. Assuming, on heuristic physical 
basis, that G$\t) G L 2 , we can apply the Titchmarsh theorem [37, 38, 39, 40] and 
deduce that the linear susceptibility x£ 0^) is a holomorphic function in the upper 
complex w-plane and the real and the imaginary part of are connected to each 




(5) 
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other by Hilbert transform. 

According to a general property of Fourier transform we know that the short term 
behavior of G$\t) determines the asymptotic properties of x& (u)- We shall obtain a 
more quantitative result by exploiting that: 

r+oo i / ■ \ j(fc+l) 

/ dte(t)t k exp[iuot] = (-i) k -j- P- + 7r5(ou) « A;!——- (6) 
J-oc du \ u J u( k+1 > 

where in the second equality we have neglected the fact that the solution is a distribu- 
tion and considered w^O. Therefore, if the Taylor expansion of the Green function in 
the limit t — > + is of the form: 

G l i\t) ^aS(t)t p + o(t p ) (7) 

the high frequency behavior of the linear susceptibility, i.e. the limit u — > oo, is: 

X^Cw) waw-^ + oCw-"- 1 ) (8) 

where a = ai^ +1 ^(3\. The parameters (3 (which is an integer number) and a depend 
on the observable $, on the specific features of the unperturbed system, and on the 
forcing under consideration. Taking into account (5) and assuming that u is real, we 
obtain that x$ — [x^\~ w )]*) so that Re[%] is an even function while Im[x] is odd 
function of u. Thus a = an is real if /3 is odd, whereas a = iai is imaginary if (3 is 
even. 

Taking into account the Titchmarsch theorem, using that X$\u) = [x$\— 
and considering the asymptotic behavior of the susceptibility, it is possible to show 
that the real and imaginary part of the linear susceptibility obey the following set of 
general KK dispersion relations [40]: 



POO 

■V^Imfx^HHP / di/ 
* Jo 



^Re^V)] (g) 



u 2 — UJ 2 



Z J Q V — UJ 

with P indicating integration in principal part and p = 0, . . . , (j3 — l)/2 if j3 is odd and 
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p = 0, . . . , /3/2 if /3 is even. Note that the faster the asymptotic decrease of the suscep- 
tibility, the higher the number of independent constraints due to KK relations it has 
to unavoidably obeys. As thoroughly discussed in [9], in the case of quasi-equilibrium 
system, the fluctuation-dissipation theorem ensures that the imaginary part of the 
susceptibility describing the response of a given observable to a perturbation is pro- 
portional to a suitably defined power spectrum in the unperturbed system. Therefore, 
observing the unperturbed system and using Eq.(10) it is possible to reconstruct the 
entire linear susceptibility, and so know everything about the response properties of 
the system. In the case of a non-equilibrium system, as discussed in the Introduction, 
this procedure is not possible. 

It is possible to use the KK relations to define specific self-consistency properties 
of the real and imaginary part of the susceptibility. We first consider the following 
application: we set p = in Eq. (10) and take the limit u> — > 0. We obtain that for 
any observable: 



which says that the static susceptibility (i.e., in a more common language, the linear 
sensitivity of the system) is related to the out-of-phase response of the system at all 
frequencies. In other terms, Eq. (11) is an exact formula for the linear susceptibility 
of the system. Note that the static susceptibility is a real number because, thanks to 
the symmetry properties discussed above, Im[%W(0)] = 0. Moreover, we know that 
Re[x^(0)] is finite because the susceptibility function is analytic (and so in particular 
non singular). This is consistent with the fact that as u — > 0, the imaginary part 
of the susceptibility goes to zero at least as fast as a linear function, as only odd 
positive integer exponents can appear in its Taylor expansion around oj = 0), so that 
the integrand in Eq. (11) is not singular. Similarly, we obtain that for u ~ the real 
part of the susceptibility is in general of the form c\ + C2U). 2 + o(w 3 ), where C\ and Ci 
are two constants and c\ is exactly given by Eq. (11). 

By exploring the u — > 00 limit in Eqs. (9)-(10) we obtain further integral con- 
straints. By applying the superconvergence theorem [54], we obtain the following set 
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of vanishing sum rules (see [39] ) : 



I 



oo 



dr/^ImfxgV)] =0 



< p < 0/2 - 1, 
0<p<(/3-3)/2, 



even 



odd 



(12) 



/ 

Jo 



oo 



V 



2p Re[ X l i\u)]diy = 



0<P< P/2-1, 
0<p<(/3-l)/2, 



/3 
/3 



even 



odd 



(13) 



Note that if (3 = no vanishing sum rules can be written for the susceptibility, whereas 
if (3 = 1 only Eq. (13) provides a zero-sum constraint. For each set of KK relations, 
an additional, non-vanishing sum rule can be obtained. If /3 is odd, the non-vanishing 
sum rule is: 



where the a constants are defined in Eq. (8). These sum rules provide additional 
general constraints that must be obeyed by any system and can be used to test the 
quality of the output of any model wishing to describe it. If we are not in the perfect 
model scenario (e.g., we use a simplified representation of some degrees of freedom) 
the sum rules can in principle be used to provide a fit for the parametrization. 

We underline that it is possible to generalize the KK theory for specific classes 
of nonlinear susceptibilities for both quasi-equilibrium and non-equilibrium systems. 
Such results, which are particularly suited for studying the fundamental properties of 
harmonic generation processes, are thoroughly discussed in [17] and will not be reported 
here. 




(14) 



while if /3 is even, we have: 




(15) 
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2.3 A practical formula for the linear susceptibility and con- 
sistency relations between susceptibilities of different ob- 
servables 

As discussed above, the definition of the linear susceptibility does not depend on the 
function f(t) modulating the additional forcing, so that it is possible to draw general 
conclusions on its properties even by choosing a specific function f(t). 

Let's consider f(t) = 2ecos(wt). The impact of the perturbation on the evolution 
of a general observable is defined as: 

5$ e (t, t , x ) = $ 6 (i, to, x ) - $ (t, t , x Q ) (16) 

where Xq and to are the initial condition and the initial time, and we associate the 
lower index e to the strength of the forcing. The Ruelle's response theory ensures that 
6<& e (t,t ,Xo) = 0(e). Following [41, 43] the linear susceptibility results to be: 

X* ( w ) = lin ^ J im xi 7(w,x ,e,T) (17) 

where: 

i i 

xi 1} (w,x ,e,T) = -- / dt<5$ £ (t,x )exp(Lvt) (18) 
1 e Jo 

is the total susceptibility, affected by the finite time and finite size response of the 
system. This quantity depends on the initial condition and in principle contains infor- 
mation about the response of the system at all order of nonlinearity. 

Since d(5$ e (t, x ))/dt = 5(d<& e (t, x )) /dt, thanks to the linearity of the time deriva- 
tive, by considering Eqs. (17)-(18) and performing an integration by parts, we obtain 
that [41]: 

x«M = -iu, x i 1} M- (19) 

Let's now find a different expression for ^ (w). The time derivative of $ in the 
unperturbed system is 

= T(x), (20) 
where T = F ■ V3>. Similarly, the time derivative for $ in the case of the perturbed 
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motion described by x = F{x) + X(x)f(t) = F(x) + 2e cos(ut)X(x) is: 



$(x) = T(x) + 2e cos(wt)H(x), (21) 

where H = X ■ V$. From Eqs. (20)-(21) we obtain that 

5$ e (t, x ) = <5r 6 (<, a*) + 2e cos(wt)S(t, x ), (22) 

where all terms are of 0(e). Furthermore, we integrate each term in Eq. (22) as in Eq. 
(18), take the limits e — > and T — > oo, and obtain: 

"T 1 1 rT 



l l r u f 

lim lim — - / di5$ e (t, xo) exp(icut) = lim lim — - / dt<5r e (t, Xo) exp(io;t) + 
e->OT->ooT e £-+or->ooT e ,/ 

1 1 ( T 

+ lim lim — — / dteS e (t, xo) [expnatf) + exp(— kjt)] expnatf). (23) 

e->0 T->oo T e Jq 

Using the definition in Eq. (17) and the identity given in Eq. (19), we derive: 

1 f T 

-ivx*\v) = Xr\u)) + lim lim — / dtS e (t,x ) + 

1 1 f T 

+ lim lim / dteEJt, x ) exp(2iwt). (24) 

e^O T->oo T e Jo 

The first limit in Eq. (24) gives, by definition, (S)o, whereas the second limit vanishes 
as the expression under integral is 0(e 2 ), since it is related to second order harmonic 
generation nonlinear process [43]. Concluding, we obtain the following general consis- 
tency relation for the linear susceptibility: 

xi\uj) = -x?(u) + -(Z)o. (25) 

UJ UJ 

Such an identity related the susceptibility of an observable $ to the susceptibility of 
the projection of its gradient along the unperturbed flow V and to the average value 
in the unperturbed state of the projection of its gradient along the perturbation flow. 
Note that the two terms on the right hand side are radically different. Whereas the 
first term is related to the projection of the dynamics along the unstable manifold, 
the second term depends on the structure of the forcing X(x), which may be entirely 
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unrelated to that of the unstable manifold. This is the fundamental reason why the 
fluctuation-dissipation theorem does not apply in the non-equilibrium case. 

Moreover, since, as shown in Eq. (7)-(8), the susceptibility of a generic observable 
decreases to zero at least as fast as w" 1 , for large values of to we have that x$ ( w ) ~ 
i/u)(E)o unless u(E) = 0. If lo(E)o ^ 0, we also have that the leading order of the 
short-time expansion of the Green function is of the form: 

4 1) (t) = 6(t)(S)o + o(t°), (26) 

in agreement with what can be found by direct inspection of Eq. (3). 



3 Application of the Response Theory to the Lorenz 
96 Model 

3.1 Statistical properties of the unperturbed Lorenz 96 Model 

The Lorenz 96 model [45, 46, 47] describes the evolution of a generic atmospheric 
variable defined in N equally spaced grid points along a latitudinal circle and provides a 
simple, unrealistic but conceptually satisfying representation of some basic atmospheric 
processes, even if such one-dimensional model it cannot be derived ab-initio from any 
dynamic equation via subsequent approximations. The evolution equations can be 
written in a scaled form as follows: 

= £i_i(> m - Xi_ 2 ) ~ Xi + F (27) 

where i = 1,2, , N, and the index i is cyclic so that X{ + n = %i-N = The 

quadratic term in the equations simulates advection, the linear one represents thermal 
or mechanical damping and the constant one is an external forcing. The evolution 
equations are invariant under i — >■ i + 1, so that the dynamics is the same for all 
variable. The time scale of the system is given by the damping time, which corresponds 
to five days. The L96 system shows different features, as different choices of F and N 
may strongly alter the topology of the attractor, alternating periodic, quasi-periodic 
and chaotic behavior in a non trivial way. However with a suitable choice of the 
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parameters N and F, the system is markedly chaotic. In particular, as F controls the 
energy input into the system, we expect that for relatively high values of this parameter 
the system should simulate a turbulent behavior and live on a strange attractor. As 
an example, setting N = 40 and F = 8 the system possesses 13 positive Lyapunov 
exponents, the largest corresponding to a doubling time of 2.1 days, while the fractal 
dimension of the attractor [21] is about 27.1 [46]. 

When computing the time derivative of the total energy of the system, defined as 
E = 1/2 ^i; the advection terms cancel. The evolution equation for E results to 
be: 

E = -2E + F^Xi. (28) 

i 

As the dynamics takes place inside a compact set, ^(xj, . . . ,xn) is bounded for any 
choice of the function \&. Therefore the ensemble mean with respect the po (o r time 
average) of the temporal derivative ^ vanishes. Therefore, defining M = Xi as the 
total momentum of the system, we obtain the following identity: 

2 (E) = ( x i) = FJ2 (*t)o = F < M >o • (29) 

i i 

Similarly, we can deduce an additional consistency relation by investigating the expres- 
sion of the time derivative of M: 

(M) = iVF+(C 2 ) -(C3) (30) 

where (C 2 ) = J2i ( x i x i-2) an d (£3)0 = J2i ( x i x i-3)o- Higher order consistence rela- 
tions can be obtain in a similar fashion. 

The equivalence of all the variables implies that over the unperturbed flow each 
observable O of the whole system satisfies ^ (O(xi)) = N {0(xj)) Q Wj. Therefore, 
we define the average energy per grid point e(A r , F) and the average momentum per 
grid point m(N, F) as: 

e 
m 

where the choice of % is arbitrary and the and F-dependence is dropped for shortness. 



2 

( x i/0 



(E)o 
N ' 

N ' 



(31) 
(32) 
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Defining q = (Ci) /N, Co = 2e, we can rewrite Eqs. (29)-(30) as follows: 



2 e = Fm, 



(33) 



m = F + c 2 - c 3 . 



(34) 



Expressing either the average energy e or the average momentum m per grid point as 
a function of the two free parameters N and F would allow to get a closure equation 
for the statistical properties of the unperturbed Lorenz 96 system. We have computed 
e(N, F) and m(N, F) by performing long integrations for values of F ranging from 
6 to 50 with step 1 and for values of N ranging from 10 to 200 with step 10. In 
all of these cases, chaotic motions are observed. We have consistently found that, 
within 0.5%, e(N,F) = e(F) and m(N,F) = m(F), so that they can be considered 
intensive quantities. Therefore, we can interpret Eqs. (33)-(34) as equations providing 
a definition of the thermodynamics of this simple one- dimensional model of atmosphere. 

The F-dependence of e and m can be closely approximated in terms of simple power 
laws. We obtain, within a precision of about 1% in the considered domain, that 



with A ~ 1.15 and 7 ~ 0.35. Such a smooth dependence of the intensive energy and 
momentum with respect to the forcing parameter F is indeed in agreement with the 
hypothesis that the invariant measure is deformed in a very regular fashion not only 
locally, but over a large range of the parameter's space. 

Note that, at the fixed point of the system corresponding to a purely zonally sym- 
metric dynamics (xj = F, Vi) we have m = F and e = F 2 /2. These formulas give much 
higher values for both m and e than what found with our empirical power laws for the 
attractor in the chaotic regime. In fact, at such an equilibrium, which is unstable in 
the parametric range explored here, the energy dissipation is much weaker than in the 
co-existing chaotic attractor, which corresponds to the case where breaking nonlinear 



m{F) = AF 7 , 



(35) 



and, consistently with Eq.(33), 



e(F) = -F 1+7 , 



(36) 
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waves and turbulent motions are present. Interestingly, the presence of well-defined 
scaling laws with respect to the forcing parameters for the energy and momentum of 
the system with different characteristic exponents in the chaotic regime and in the co- 
existing unstable equilibrium is in agreement with previous finding recently obtained 
in a simple baroclinic quasi-geostrophic model [53]. 

3.2 Asymptotic properties of the linear susceptibility 

We perturb the L96 model by adding a small perturbation modulated by f(t) = 
2ecos(W). The resulting evolution equation is: 

d x ' 

-jy = - x i+2 ) - Xi + F + 2e cos(ut)Xi (37) 

where Xi = Xi(xi, . . . , x^) is a generic function of variables X{. We adopt the chaotic 
hypothesis [23, 22] and we follow the theory proposed by Ruelle [10] and discussed in 
Sect. 2 in order to study the linear response of suitably defined observables to the 
perturbation. We first propose to study the high-frequency, response by analyzing in 
detail the asymptotic properties of the resulting susceptibilities. As discussed in section 
2, this constitutes a crucial step for constructing the set of applicable KK relations and 
for computing the value of the sum rules. 

We consider two different forcing patters Xi. In the first case, we apply the pertur- 
bation over all the grid points and we choose Xi = IV %. Given an observable $, we 
refer to the linear Green function and linear susceptibility resulting from this choice 
of Xi as G$ N and X* jv> respectively, where the lower index N indicates that the per- 
turbation acts over all the variables x%. We refer to this pattern of forcing as global 
perturbation. In the second case, we apply the perturbation only on the variable Xj 
of the L96 model, and we choose X 1 = OVz ^ j, Xj — 1, . Since all the points are 
equivalent in the unperturbed case, the choice of j is arbitrary. In this case, when 
referring to the linear Green function and the linear susceptibility, the lower index 1 
substitutes the N, indicating that the perturbation is localized to one point. We refer 
to this pattern of forcing as local perturbation. 
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3.2.1 Global perturbation 

We consider perturbations with spatial pattern given by Xi = 1 V i and analyze the 
response of the observable E. Following Eq. (3), the linear Green function Gg a (t) can 
be explicitly written as: 



G%{t) = I p Q (dx)e(t)AU (t)E(x) 
Po (dx)e(t)l-VE(x(t)) 
J Po (dx)Q(t)J2^(E(x(t))) (3* 



where x(t) satisfies the unperturbed evolution Equation (27). Taking into account Eq. 
(7) and Dq. (8), in order to obtain the asymptotic behavior of the susceptibility, we 
need to study the short time behavior of the Green function. Therefore, we express 
E(x{t)) as a Taylor series about t = considering the unperturbed flow, compute 
the integral of each coefficient of the t-expansion over p , and seek the lowest order 
non- vanishing term [43]. The first two terms of the Taylor expansion of E in Eq. (38) 
give: 



G { i] N (t) = J p o (dx)0(t) J2 d i (^l*=o + tE\ t=Q + o(t) 

= J Po (dx)e(t)[J2 x i- {J2 2x *- NF ) t + °^ 



(39) 



Using Eqs. (7)-(8), the leading terms of the asymptotic behavior of linear susceptibility 
can be written as: 

xSvH = * ( E fc>o ) / w + ( E < 2 ^>o - NF ) i" 2 + °(^" 2 ) 

i i 

= iN N — + cj" 2 40 

Since the symmetry with respect the index % is valid also in the perturbed case, given 
our choice of the forcing pattern, the linear susceptibility of the total energy is given 
the sum of N identical contributions, each corresponding to the susceptibility of the 
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observable e = l/2x? for each of the N variables Xi of the system. Therefore, it 
is possible to define an intensive linear susceptibility x^n = ^-/^Xeni wriere XeN 
describes the response of the local energy to the external perturbation. In particular 
in the limit w->oowe have: 

X ( eM U ) = % — + °( W )• ( 41 ) 

UJ UJ 

Equations (40) and (41) imply that the imaginary part dominates the asymptotic 
behavior of the susceptibility, so that at high frequency the response is shifted by 
about 7i / 2 with respect the forcing. Observing that the leading term of asymptotic 
X is of order uj' 1 just one sum rules apply for either susceptibilities. Limiting our 
attention to the intensive quantity e, by applying Eq. (15) we obtain: 

/CO 
Re[ X «M] = - m. (42) 

Along similar lines, if we select as observable the total momentum M, we derive that 
the asymptotic behavior of its linear susceptibility is: 

X { mU") = Nx%(f») =i~ + - 2 + o(^" 2 ), (43) 

where we have defined the intensive susceptibility X^iv^)' where \x is the intensive 
momentum of the system. As in Eq. (41), the asymptotic behavior is determined by 
the imaginary part of X-, an d the real part of the susceptibility provides the following 
sum rule: 



I 

Jo 



Re[xSr(w)]dw = ;. (44) 



7T 



We now wish to go back to the general consistency equation for linear susceptibilities 
given in Eq. (46). Considering that in the perturbed system the time derivative of the 
total energy of the system can be written as: 

E = -2E + FM + 2e cos(ujt)M (45) 
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the general result given in Eq. (25) can be written as follows: 



since in this case T = —2E + FM and H = M. It is easy to check that the asymptotic 
behavior for the susceptibilities given in Eqs. (40)-(43) is in agreement with Eq. (46), 
which is valid at all frequencies. 

3.2.2 Local perturbation 

We now perturb the system in a single grid point. The symmetry of the unperturbed 
system implies the equivalence of every point of the latitude circle. Indicating with Xj 
the grid point where forcing is exerted, the pattern of the perturbation vector field is 
Xi = OVz 7^ j, Xj = 1. We consider the same modulating monochromatic function 
f(t) = 2ecos(ut) as in the previous case. We analyze the asymptotic behavior of the 
linear susceptibility of the total energy of the system E, of the total momentum of 
the system M, which are global variables, and of the local variables constituted by the 
energy Ej = l/2x| and momentum Mj = Xj of the perturbed grid point and of its 
immediate neighbors. 

Since we are looking at the linear response and the global perturbation is given by 
N spatially shifted copies of the local perturbation, for any observable $ of the form 
...,x n ) = 4>(xi) we must have xji = X*/ N jf = An- 

In the case of the observable E, it is straightforward to verify the previous identity 
at least in the asymptotic regimes. In fact, the short time behavior of the Green 
function describing the response of the E to the local perturbation results to be: 

Ge,M = J Po(dx)e(t)d j (E\ t=0 + tE\ t=0 + o(t)) 

= j Po (dx)e(t)d,[\j2^- (J2 2x i-^ F y_ 

= 0(t) ( (2x,) - (2 Xj -F) t + (t)) , (47) 

so that the asymptotic behavior of the corresponding linear susceptibility susceptibility 
is: 

X%i = i — + o(^ 2 ), (48) 
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which agrees with what found for the intensive energy response when the global per- 
turbation is applied (see Eq. (41)). The sum rule for the real part of the susceptibility 
is exactly the same as in what given in (42): 



Re[xS(u,)]dc = -. 



(49) 



Analogously, we obtain that the asymptotic behavior of x$ e i can ^ e wr hten as: 



*2» = i- + ^ (50) 

UJ UJ 

with the corresponding sum rule: 

He|xff»]d W =- ) (51) 

in perfect agreement with Eqs. (43) and (44), respectively. 

It is rather interesting to look into local energy observables. Considering the energy 
Ej of the perturbed grid point Xj we have that its short term Green function can be 
written as: 

G^iW = J Po(dx)Q(t)dj ^xj + (xjixj^Xj+i - x i _i^_ 2 - Xj + F)^jt + o(t) 

= ©(*) (xj) - (xj-iXj+i - Xj-iXj-2 - xj + F) t + o(t) . (52) 

Since 

= (^)o = (xj-ix j+1 - Xj-xXj-2 ~ Xj + F) Q (53) 

because the po-average of the temporal derivative of any observable vanishes, thanks 
to the compactness of the attractor, we obtain the following asymptotic behavior for 
the linear susceptibility 

XeU=*- + -2+°^ )■ ( 54 ) 
J UJ UJ 

Since, by linearity, Xe\ = EaXe n comparing this result with what obtained in 
Eq. (48), we note that the susceptibility of the energy at the position of the forcing 
Ej provides the leading asymptotic term to the susceptibility of the total energy E. 
Consequently, in the high-frequency range Xe- i ~ Xeii an< ^ the two susceptibilities 



1 1 
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obey the same non-vanishing sum rule, so that: 



/oo 
Re[ X g>)]du; = -m (55) 

Nevertheless, by comparing Eqs. (48)-(54), we discover that contributions to the second 
leading order (oc oj~ 2 ) in the high frequency range of the susceptibility of the total 
energy do not come just from the response of the energy at the perturbed grid point, 
perturbed grid point but some other point give a contribution of order uj~ 2 . Therefore, 
the asymptotic behavior of the real part of Xei * s no ^ captured by Xe v The locality 
of the interaction suggests to look at the energy of the closest neighbors of Xj. Because 
of the asymmetry of the nonlinear terms in the L96 evolution equations, we consider 
the observable ip = l/2(Ej +1 + Ej +2 + Ej_i). It is possible to prove that: 

4» = - <*" (J *' - ^ = -<^> + o(^). (56) 

It is easy to observe that the sum of x^\ an d Xe] i P rov ides the correct leading order to 
the asymptotic behavior of both the real and imaginary parts of Xei- We shall provide 
an argument why this strongly supports the close resemblance of the two functions 
X { e\ and XeI c ,v where E ioc = Ej + ij) = \j2x) + l/2x| +1 + l/2xf +2 + l/2x J 2 _ 1 is the 
energy of the cluster of points centered in xj. 

The analysis of the asymptotic behavior of the susceptibilities related to the local 
momentum of the system provides additional insights. It is possible to prove that for 
large frequency the linear susceptibility of the momentum of the perturbed grid point 
is: 

^>) = i- + A + o(uT 2 ), (57) 



which suggests that the response of the local momentum captures the correct asymp- 
totic behavior of both the real and the imaginary part of the total momentum M. 
Concluding, we obtain that the following sum rule can be stablished: 



oo 

Re[ X «>)]du; = -. (5? 



Therefore, such a constraint is exactly the same whether we analyze the the response of 
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the momentum of a single variable when the perturbation acts over all the grid points, 
or of the total momentum in the case of a local perturbation, or, in this latter case, of 
the momentum of the grid point where the local perturbation is applied. 

As we have seen in this section, the coefficients of the leading asymptotic terms 
and the sum rules are given by simple linear functions of m (or equivalently, thanks 
to Eq. (33), by e) and by F. As we have proposed an efficient parameterization of m 
and e as functions of F alone in subsection 3.1, our results can be easily applied and 
numerically verified for a very large class of L96 models. 

4 Results 

4.1 Simulations and Data Processing 

The accurate calculation of the linear susceptibility of the general observable $ is not 
as easy task, since the definition given in Eq. (17) requires the evaluation of two limits, 
whereas we can actually compute only the quantity given in Eq. (18). Averaging the 
response over a long time T allows for improving the signal-to-noise ratio. Noise is 
present because, due to the chaotic nature of the flow, we have a continuous spectral 
background. Instead, considering small values for the perturbation strength e degrades 
the signal-to- noise ration, but, on the other hand, it is crucial to select a small e in 
order to keep the perturbations as close as possible to the linear regime. As discussed 
in [43], we can improve the signal-to- noise ratio without needing to perform very long 
integrations and to consider large values for e by performing an ergodic averaging of 
the quantity averaging the quantity x^i 00 ^ x h e > T). Therefore, we choose the best 
estimator of the true susceptibility x$ as: 



where the Xi are randomly selected initial conditions chosen on the attractor of the 
unperturbed system. 

The numerical integrations of the Lorenz 96 system have been performed using the 
standard configuration where N, the number of degrees of freedom, is set to 40, and 




(59) 



i=l 
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F, the intensity of the unperturbed forcing, is set to 8 [47, 45]. Equations (27) and 
(37) are solved using the standard fourth order Runge-Kutta numerical scheme. 

For a given observable $, the susceptibility at angular frequency u is computed by 
applying Eq. (59) to K outputs of Eq. (37), each starting with a different initial condi- 
tion, where the perturbation has the same angular frequency u. The angular frequency 
oj ranges from ui = 0.2n to Uh = 207r with steps of 0.0l7r. Each simulation performed 
with a perturbation of angular frequency u runs from t = 0uptot = T = 400tt/uj, 
which corresponds to 200 full periods of the forcing. The length of the simulations 
depends on the corresponding period of the forcing because we are interested in ob- 
taining a frequency-independent quality for the signal. We have observed that the 
linear response approximation is obeyed to a good degree of approximation for up to 
e « 1, which implies that the third order nonlinear effects are relatively small. See 
[17, 43] for further clarifications on this latter point. 

When considering the susceptibilities describing the response to the global pertur- 
bation, we present results obtained using e = 0.25 and averaging over K = 100 random 
initial conditions. When assessing the linear response to the local perturbation, a 
reasonably clear signal is obtained using e = 1 and averaging over K = 300 initial 
conditions. 

Note that, since we are interested in the linear response, it is could have been pos- 
sible to compute the susceptibility using a generic modulating function f(t) (see Eq. 
(4)) rather than having to resort to multiple monochromatic perturbations. Neverthe- 
less, for reasons of clarity, and for emphasizing that chaotic dynamical systems can 
be analyzed using tools typical of spectroscopy, we have used a more cumbersome but 
probably more convincing approach. 

We underline that the numerical results have been obtained using a commercial 
laptop rather than resorting to HPC. This comes from the motivation of showing that 
the methodology presented is robust enough that relatively low-end means allow us 
to see the physical and mathematical properties of our interest. We emphasize that, 
using HPC, it is rather easy to greatly increase the quality of the signal by increasing 
K and/or T by a one or two orders of magnitude. 
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4.2 Global Perturbation 

We first consider x^n = V-^Xsiv> where e — E/N, and follow up from the discussion 
in subsection 3.2.1. The measured real and imaginary parts of the susceptibility are 
depicted with the black lines in Fig. (l)a,b. The imaginary part has a broad spectral 
feature (with two distinct peaks) spanning from u « 2 to u ~ 4, which corresponds 
to about twice the time scale (=1) of the system and to four times (see Eq. (28)) 
the relaxation time of the energy. This hints at the fact that it is not obvious to 
constrain the spectral features of the response an observable just by performing a scale 
analysis of its evolution equation. For higher values of u, the imaginary part decreases 
in a very regular way, so that in the upper range a very good agreement with the 
asymptotic behavior ~ m/u presented in Eq. (8) is obtained. For low frequencies, the 
imaginary part appears to decrease towards zero, as expected from symmetry reasons. 
Instead, the real part presents a dispersive structure in correspondence with the broad 
maximum of the imaginary part, and changes sign for u ~ 6, so that it is negative for 
high values of the frequency range. The asymptotic decrease to zero in this range is also 
in excellent agreement with the estimate ~ — (F — 2m) /u 2 given in Eq. (8), whereas 
for low frequencies the real susceptibility tends to a very high value, this suggesting 
that the strongest response is obtained for static perturbations. 

The measured real and imaginary parts of X^jv = ^-/^Xmnj w here /i = M/N, are 
depicted in black in Fig. (2)a,b. Interestingly, the spectral feature of the imaginary 
part is shifted to higher frequencies than in the case of the energy susceptibility, so 
that a well-distinct peak centered on value of u ~ 6, which approximately corresponds 
to the natural time scale of the system. For low frequencies, the susceptibility has 
almost exclusively a real component. As opposed to the previous case, the largest 
value for the in-phase response is not obtained for ultralow frequencies, but rather for 
u ~ 4. The asymptotic behavior of both the real and imaginary parts is in perfect 
agreement with the theoretical result given in Eq. (32), so that they are found to 
decrease asymptotically for high frequencies as 1/u 2 and 1/oj, respectively. 

We apply the truncated KK relations to the measured data to test the quality of 
the data inversion process. The estimates of the imaginary part (starting from the 
measured data of the real part) and of the real part (starting from the measured data 
of the imaginary part) obtained by applying Eqs. (9)-(10) are shown for XeN m ^ ue 
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in Fig. (l)a,b and for x\tN m Fig- (2)a,b. We observe that whereas agreement is very 
good for the real part for both susceptibilities, only a qualitative match is obtained for 
the imaginary part, with large discrepancies for uj < 2. In this latter case, moreover, 
the well-known problem of KK divergence at the boundaries of integration [39, 40] is 
very serious for uj = UJi. 

It is crucial to test whether the discrepancies are due to the fmiteness of the spectral 
range or are, instead, due to basic problems in the applicability of the Ruelle response 
theory, related to the fact that the invariant probability measure of the unperturbed 
system actually features large deviations from an SRB measure. 

We proceed testing the first case. In order to widen the spectral range over which the 
susceptibility is defined, we will exploit the asymptotic properties obtained in section 
3.2 as well as the low frequency behavior of the susceptibility discussed in section 2. 
We redefine the the imaginary part of the susceptibility of Xen as follows 

f %Im[x { S{ui% 0<u<u t , 
lm[ X ^ N (oj)\ = I Mx^u;)], uj l < uj < uj h , (60) 

w > w fcl 

where the measured data are sandwiched between the low and high frequency limit. 
Whereas we have a rigorous result for the high frequency limit, the low frequency 
limit is computed by making the reasonable assumption that the leading order of the 
uj — y limit is linear (see discussion after Eq. (11)). Similarly, the real part of the 
susceptibility x^n can be redefined as follows: 

< uj < oji, 

ui <uj < u h , (61) 
uj > UJ h , 

where we have used the fact that at low frequencies the real part of the susceptibility 
is constant in w up to a quadratic term. A corresponding procedure is used to extend 
the spectral range of the X^n( u ) > where the suitable asymptotic behaviors described 
in subsection 3.2.1 are adopted. The red lines in Figs. (l)a,b-(2a,b) present the results 
of such extrapolations, and the magenta lines show the outcome of the data inver- 
sion of these functions performed via KK relations. We observe that the agreement is 



MA 



UJ)] 
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outstanding, with almost perfect overlap inside the region where measurement is per- 
formed and remarkable agreement also in the low and high frequency range. This is a 
very convincing evidence that the Ruelle response theory can be successfully applied for 
this system. Since the KK relations provide, first and foremost, consistency tests, the 
agreement the original and the KK-transformed susceptibility automatically confirms 
that the extrapolation procedure we have adopted is correct. A still better agreement 
would be found had we taken into account value of u larger than what considered in 
the extrapolation used here (up to 100 n). 

Furthermore, let's consider the results presented in subsection 3.1. The slopes of 
the functions e(F) and m(F) are given by 



dF 



A7F 7 , (62) 



!£(£) = xH + tIf, = m(F) (ill) . (63 ) 

dF 2 V ' 2 K ' 

They correspond, by definition, to the static susceptibility of the observables e and m, 
respectively, for the global perturbation with = 1 considered here. When evaluating 
the derivatives of e(F) and m(F) for F = 8 we obtain (d e(F) / d F) p =8 ps 1.6 and 
(dm(F)/dF)F=8 ~ 0.11. These values are in good agreement with what found by 
extrapolating the corresponding real part of the susceptibilities for u — > via KK 
relations and shown in Figs, (la) and (2a). 

Apart from the verification of the validity of KK relations, we want to provide 
further support for the quality of the linear susceptibilities considered. 

First, we test the sum rules given in Eq. (42) and (44) for the real part of the 
extrapolated susceptibilities xflii 00 ) an d xL^ivO^); respectively. Our findings are pre- 
sented in Fig. (3, where it is shown that an excellent agreement (within 1%) is found 
between the theoretical values and the numerical results. Since Re[xe*e,a](w) is negative 
in the high-frequency range, the convergence of the integral to the theoretical value of 
the sum rule is from above, whereas the opposite occurs for Re[ X $ eia (uj)]. Extending 
the integral for even larger values of u would bring the numerical results to an almost 
perfect agreement with the theory. 

Following the definition given in Eq. (3), the Green function G$\t) computed for 
an observable $ and a given pattern of perturbation flow -Xj(x) (in this case Xi = lVz) 
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can be used to compute the time- dependent linearized impact of all perturbations 
with the same spatial pattern Xi(x) but with arbitrary time modulation. Whereas 
the direct estimate of the Green function from the time dependent dynamics can be 
obtained by performing an ensemble of simulations where the time modulation of the 
perturbation is given by a S(t) pattern (see discussion in section 2, we take the indirect 
route by considering Eq. (5). By applying the inverse Fourier Transform, we derive the 
Green functions corresponding to Xe]l(uj) and x^a(uj). The results are presented in Fig. 
(4): for both observables the Green functions are clearly causal, and their short-time 
behavior agrees remarkably well with what be deduced by looking at the asymptotic 
properties of the corresponding susceptibilities (compare Eqs. ((41) and (43)). 

4.3 Local Perturbation 

The data obtained for the numerical simulations of the response to the local perturba- 
tion are, given the much weaker overall strength of the forcing, much noisier that those 
presented in the previous section. Nevertheless, we shall see that all the theoretical 
predictions are verified to a surprisingly good degree of approximation. 

The global observables E and M are of the form <3?(xi, . . . , x n ) = Yli=i *&( x i)i where 
<f)(x) = x 2 /2 and <&(x) = x, respectively. We have consistently verified that the identity 
= X$/jvjv = X$,N discussed in the previous section applies in the whole spectral 
range explored by our simulations, compatibly with the (slightly) different signal-to- 
noise ratios in the two sets of simulations. See, e.g., Fig. 5 for the comparison between 
the two susceptibilities Xei anc ^ Xsn- 

We then proceed to analyze more in detail the linear susceptibilities related to local 
observables. In Fig. (6) we present our results concerning the real and imaginary part 
of x^EjA 00 )- Analogously to what observed in the previous subsection, we have that 
once the measured susceptibility is extrapolated using the theoretical results obtained 
via response theory and KK relations, we have an excellent agreement between the 
original real and imaginary parts and those obtained using the KK inversion. The 
KK algorithm, instead, provides only a partially satisfying outcome when only data 
from the measured range are considered. Relatively discrepancies are found near the 
boundaries of the data range, with an especially serious bias near ui for the imaginary 
part of the susceptibility. 
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When comparing Xe],!^) anc ^ Xe\( u ) ( see Fig- 5), we observe that for low fre- 
quency the response of the energy at the grid point where the perturbation is applied 
accounts for about half of the response of the total energy, thus implying that the re- 
maining half is redistributed among the remaining N — 1 grid points. The relevance of 
the response of grid points other than the directly perturbed one also explains why the 
peak of IjxiXe\((jj) (and so of lmx^ N (w)) than that of Imx^-^aj) - see the frequency 
range 2 < u < 4. Slower perturbations allow other grid points Xk ^ Xj to respond 
effectively. 

Instead, since the leading asymptotic order of Xe an< ^ Xei( u ) * s ^ ne same > a t 
high frequencies the local energy response accounts for most of the energy response of 
the whole system. In this case, the incoming perturbation is so fast that the internal 
time scales of the system as bypassed, and mainly a local effect is observed. Neverthe- 
less, the second leading order of the asymptotic expansion of the two susceptibilities 
has opposite sign (see Eqs. (41) and (54)), which suggests that at any large but finite 
frequency the local energy response is only a good approximation to the response of 
the total energy. The changeover between the two regimes occurs around the frequency 
of the peak of Imx^ 1 (w), which corresponds to a perturbation with period close to 1. 

Thanks to the asymptotic equivalence between Xe] an d Xe\(°°) ( an< ^ X^iU^))' 
they must obey the same sum for the real part of the susceptibility (see Eqs. (42)- 
(55)), even if the two real parts, as discussed above, are rather different in value in the 
low frequency range and even in sign in the high frequency range. Figure 8 confirms 
that this rather counter-intuitive behavior is actually observed. Note also that sum 
rules, resulting from an integration, are less sensitive to noise in the data, but this 
occurs if and only if the underlying signal is correct. Therefore, we understand that in 
X^i(w) and XeN^) t ne strong static and quasi-static response and the (rather odd) 
negative sign for high frequencies of the real part of the linear susceptibility, which are 
crucially related to the behavior for the grid points different from the perturbed one) 
compensate each other to guarantee agreement with the sum rule obtained from the 
real part of x^e- i( w )' which instead has a smaller range and more regular (monotonic) 
behavior with frequency. 

A formally similar - and analogously spectacular - spectral compensation has been 
observed in a physical process as different from what we are analyzing here as the 
electromagnetically induced transparency [55]. The result obtained here supports pre- 
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vious findings obtained on quasi-equilibrium systems suggesting that sum rules do not 
depend on many-particle interactions [39, 40]. 

The investigation of the linear susceptibility of the variable Xj is not as insightful 
as that of Ej. We find that linear susceptibility X^c is quite similar to Xmi( u ) 
(and xi atO^)) see Fig- 2) in both the real and imaginary parts at all frequency. The 
only notable differences are that the static response Xj is slightly larger than than of 
M, and that the imaginary features a secondary peak at slightly larger frequencies 
than the main spectral feature. We have verified, as in the previous cases, the results of 
the numerical simulations accurately agree with the theoretical results regarding the 
asymptotic behavior of both the real and imaginary part and that KK relations map 
to high degree of precision the real and the imaginary parts into each other. See Fig. 
7 for details. 

We present as main finding of the analysis of the observable Xj that, as predicted 
by the theory, the real part of xi^iX 6 ^) obeys the same sum rule as the real part of 
X^mi(°°) 01 °f X^ivO^)) because the corresponding imaginary parts feature the same 
asymptotic behavior. Figure 8 shows that in the case of the momentum variables the 
cumulative integral is rather similar for the susceptibility of the local and of the global 
variable, with small discrepancies in the region around the peak of the response. 

4.4 Further implications of Kramers-Kronig relations and sum 
rules 

We now show how the knowledge of the asymptotic behavior of the real and imaginary 
part and the knowledge of the validity of the KK relations and related sum rules allow 
to draw general conclusions on the similarities and differences between two given linear 
susceptibility functions. Let's consider the case that these two susceptibilities feature 
the same first order asymptotic expansion in the high frequency limit. Let's assume 
that it is an odd power of u, so that the real part is negligible for high frequencies. 
Therefore, the two susceptibilities will obey the same sum rule for, e.g. the th moment 
of the real part. 

If they agree also in the asymptotic behavior of the real part, they cannot feature 
large discrepancies in the low frequency range of the real part of the susceptibility either, 
or otherwise the agreement of the sum rules would be broken. Therefore, the real part 
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of the two susceptibilities are similar, and, as a consequence of the KK relations, the 
two imaginary parts will also be similar. 

If, instead, there is a discrepancy in the asymptotic behavior of the real part of the 
two susceptibilities, the two real parts will necessarily be rather different in the low 
frequency range, again in order to comply with the sum rule constraint. As the two 
real parts are different, the imaginary part of the two susceptibility will also be rather 
different, except, from hypothesis, in the high-frequency range. 

The first scenario envisioned here pertains to the pair of linear susceptibilities 
Xixj,i( u ) an d Xm,i( w )i whereas the second scenario is related to the pair of linear suscep- 
tibilities Xe*. i( w ) an d Xe\( u )- Note that, taking into account the asymptotic properties 
of the susceptibility of the observable Ei oc = l/2x? + l/2x^ +l + l/2Xj +2 + l/2x? 1 , dis- 
cussed in the previous section we conclude that Xe i an< ^ Xei should be similar for 
all values of u. 

Obviously, a similar argument applies if the leading order is real. This discussion 
further clarifies that the higher the number of independent KK relations and related 
sum rules verified by a susceptibility functions, the more stringent are the constraints 
on its properties. 

4.5 Additional Properties of the Linear Susceptibility 

The special mathematical properties of the linear susceptibilities allow to investigate 
further properties of the response. In particular, we note that for m > 1 the function 
[x$''] m is analytic in the upper complex w-plane just . This allows, as discussed 

in [40] to derive, in principle, an infinite set of integral relations (KK and sum rules) 
deriving just from the holomorphic proprieties of the susceptibility. As an example, we 
have considered the square of the linear susceptibility [x^jvO^)] 2 - From Eq. (41), it is 
easy to prove that the following asymptotic expansion holds for large values of u: 

bSMP = -^^ + *1. (64) 

As shown in the first panel of Fig. 9, KK relations are found to connect up to a high 
degree of approximation the real and imaginary part of [x^v^^)] 2 - Moreover, thanks 
to the asymptotic behavior given in Eq. (64), it is possible to establish, thanks to Eqs. 
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(13)-(14), the following sum rules: 



/ oVRe[ X ^H] 2 = 0, (65) 
Jo 

&vvlm[x { ^ N (v)Y = -m 2 , (66) 

The second panel of Fig. 9 shows that the obtained numerical results are in excellent 
agreement with the theoretical predictions. Note that these results do not have an 
obvious physical interpretation, as the inverse Fourier Transform of [x^jv^)] 2 is given 
by the convolution product of the Green function Gf N (t) with itself, while they depend 
only on the formal properties of the linear susceptibility. 



5 Practical Implications for Climate Change Stud- 
ies 

In this paper we have constructed and verified to a high degree of accuracy the lin- 
ear response theory for a simple yet prototypical climate model by computing the 
frequency-dependent susceptibilities of several relevant observables related to localized 
and global patterns of forcings. These results pave the way for devising a rigorous 
methodology to be used by climate models of any degree of complexity for studying 
climate change at, in principle, all time scales using only a very limited set of experi- 
ments, and for exploiting effectively the currently adopted ensemble runs methods. 

Let's consider, for sake of simplicity, that the observable $ is the time-dependent 
globally averaged surface temperature of the planet T$, that F(x) represents the whole 
set of climate equations in a baseline scenario (e.g., with pre-industrial CO2 concen- 
tration), and that the perturbation field X(x) is nothing but a constant field of CO2 
concentration, which directly impacts only the radiative part of the code. The per- 
turbation is modulated by a time-dependent function f(t) to be specified below. We 
assume, for simplicity, that the model does not feature daily or seasonal variations in 
the radiative input at the top of the atmosphere. 
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From Eq. (2), we have that 



(Tsf ] (t) = j_ 



— oo 



&tiG$ (<Ji)f(t-ai). 



(67) 



In practical terms, the left hand side of this equation is nothing but the ensemble 
average of the time series of the change between the globally averaged surface temper- 
ature of the planet at a time t after the perturbation has started. Note that the direct 
estimate of G^( a ) * s likely to be overwhelmingly difficult. Using Eq. (4), we have 
that: 



which implies that once we compute the Fourier Transform of the time series mentioned 
above and we know the modulating function f(t) (and so its Fourier Transform f{oj)), 
we can reconstruct Xts( u )- Let's select a particularly simple example of modulating 
function f{t) = e(6(t) — 6(t — r)). This is just a rectangular function of width r, of 
height e, and shifted from the origin by a forward time translation r/2. In practical 
terms, this corresponds to changing abruptly the field CO2 concentration by e at time 
t = and taking it back to its original value at t = r. we then obtain: 



we know everything about the response of the system at all time scales, including the 
static response. Note that any choice of f(t) is equally valid to set up this procedure 
as long as f(t) is square integrable. This implies that, in a very profound way, the kind 
of forcing scenarios used in the various assessment Reports of IPCC, where the CO2 
concentration typically stabilizes at a different value from the preindustrial one (so 
that f(t) does not tend to as t — > 00) are not necessarily the only nor the best ones, 
in spite of what could be intuitively guessed, to study even the steady state response 
of the system. 

Obviously, a similar set of experiments could be devised for studying rather thor- 
oughly the response of the climate system to a variety of forcings, such as changes in 
the O3 concentration, aerosols, solar radiance, as well as to changes in the parameteri- 



(Ts) (1) («)=xg(w)/(w), 





(69) 



Once we know x^j(o;), as widely discussed in this paper, we can compute G^(t), and 
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zations. In the case of uncoupled models of one subdomain of the climate system (e.g. 
atmospheric and oceanic GCMs, land-surface models), this strategy could be used to 
study the impact of perturbations to the boundary conditions provided by the other 
subdomains of the climate system. 

6 Summary, Discussion and Conclusions 

The climate can be seen as a complex, non-equilibrium system, which generates entropy 
by irreversible processes, transforms moist static energy into mechanical energy [1, 2] 
as if it were a heat engine [3, 4], and, when the external and internal parameters 
have fixed values, achieves a steady state by balancing the input and output of energy 
and entropy with the surrounding environment [56, 4]. For such basic reasons, the 
tool of equilibrium and quasi-equilibrium statistical mechanics cannot provide suitable 
tools for studying the fundamental properties of the climate system. In particular, the 
fluctuation-dissipation theorem, which allows for deriving the properties of the response 
of the system to external perturbations from the observations of its internal variability 
cannot be applied. 

It is reasonable to ask whether is possible to evaluate how far from equilibrium the 
climate system actually is. It is possible to evaluate such distance in a mathematically 
sound way by assessing the ratio of the dimensionality of the attractor of the system 
over the total number of degrees of freedom. Whereas a ratio close to one indicates that 
only small deviations from equilibrium are present, a small ratio suggests that strongly 
non-equilibrium conditions are established. See [57] for a detailed treatment of this 
problem in the classical case of heat conduction. In the case of a quasi-geostrophic 
atmospheric model forced by Earth-like boundary conditions, the dimensionality of 
the attractor of the model is about one order of magnitude smaller than the total 
number of degrees of freedom [58]. While not conclusive, this seems to suggest that 
the best framework to interpret the climate is that of a far from equilibrium system. 

Following either explicitly or implicitly the programme of the Catastrophe theory 
[59], many authors have approached the problem of understanding the fundamental 
properties of the climate system by looking at the detailed structure of the bifur- 
cations of the deterministic dynamical system constructed heuristically in order to 
represent the dynamics of the main climate modes using as few degrees of freedom as 
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possible. Such an approach often hardly allows to efficiently represent the fluctuations 
and the statistical properties of the system. The introduction of stochastic forcing 
provides a relatively simple but conceptually rich partial solution to some of these 
draw-backs, even if the Hasselmann programme [60] suffers from the need for a - usu- 
ally beyond reach - closure theory for the properties of noise. Therefore, the stochastic 
component is usually introduced ad hoc, with the ensuing lack of universality and/or 
robustness when various levels of truncations are considered. These strategies have 
anyway brought to outstanding scientific results and has been suggested the existence 
of generic mathematical structures present in hierarchies of CMs [61]. Recently, the 
unified treatment of chaotic and stochastic dynamics using the results of the mathe- 
matical theory of random dynamical systems is emerging as new, promising paradigm 
for the investigation of the structural properties of the climate system [62]. 

We have proposed a different perspective. In agreement with the view given above, 
we have taken as mathematical framework for the analysis of the climate system that of 
non-equilibrium statistical mechanics, and have focused on the steady state properties 
of ergodic dynamical systems [20] possessing the special property of having an invariant 
measure of the SRB type [21]. As proposed by the chaotic hypothesis [23, 22], this 
mathematical framework is well suited for analyzing general non-equilibrium physical 
systems. 

In this context, the impact on the system of general perturbation can be treated 
using the response theory recently introduced by Ruelle [10, 25, 16], which allows to 
compute the change in the expectation value of a generic observable as a perturbative 
series where each term is given by the average over the unperturbed invariant measure 
of a function of the phase space which depends on the considered observable and on 
the applied perturbation. In other terms, even if the internal dynamics of the system 
is nonlinear and chaotic, the leading order of the response is in general linear with the 
strength of the added perturbation. This approach overcomes the difficulties related 
to the singularity of the invariant measure discussed in [63] . 

At each order, the propagator of the perturbation, i.e. the Green function, is causal. 
This allows for applying dispersion theory and establish general integral constraints - 
KK relations - connecting the real and imaginary parts of the susceptibility, i.e. the 
Fourier Transform of the Green function [17, 43] 

In this paper we have first recapitulated the main aspects of the general response 
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theory and have propose some new general results, which boil down to consistency rela- 
tions between the linear susceptibilities of different observables. The obtained equation 
provides the basic idea why the fluctuation-dissipation theorem does not apply in non- 
equilibrium cases. 

We have showed for the first time that the Ruelle linear response theory can be 
applied with great success to analyze the climatic response to general forcings. We have 
chosen as test bed the Lorenz 96 model [45, 46, 47], which, in spite of its simplicity, has 
a well-recognized prototypical value as it is a spatially extended one-dimensional model 
and presents the basic ingredients, such as dissipation, advection and the presence of 
an external forcing, of the actual atmosphere. Such a model features a different level 
of complexity with respect to those adopted in previous numerical investigations of 
Ruelle's theory [41, 42, 43] 

We have analyzed the frequency dependence of the response of the local and global 
energy and momentum of the system to perturbations having a global spatial pattern 
and to perturbations acting only on one grid point. We have derived analytically 
several properties of the corresponding susceptibilities, such as asymptotic behavior, 
validity of KK relations, and sum rules. We have shown that all the coefficients of 
the leading asymptotic expansions as well as the integral constraints can be written as 
linear functions of parameters that describe unperturbed properties of the system, and 
in particular its average energy and average momentum. The theory has been used 
to explain differences in the response of local and global observables, in defining the 
intensive properties of the system and in generalizing the concept of climate sensitivity 
to all time scales. 

We have then verified the theoretical predictions from the outputs of the simulations 
up to a high degree of precision, even if we have used rather modest computational 
resources (a total of about 30 cpu days of a mid-range commercial laptop). We have 
verified that the linear response theory holds for perturbations of intensity accounting 
to up to about 10% of the unperturbed forcing terms. Even when local perturbation 
and local observables are considered it is possible to achieve a signal-to-noise ratio 
which permits rather satisfactory comparisons with the theory. We have proved that 
the combined use of KK relations and the knowledge of the asymptotic behavior of the 
susceptibilities allows for extrapolating in a rigorous way the observed data. We also 
have shown how to reconstruct the linear Green function, which can be used to map 
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perturbations of general time modulation into changes in the expectation value of the 
considered observable for finite as well as infinite time. 

Our numerical experiments have been performed using one of the standard settings 
of the Lorenz 96 model, namely the version identified by having N = 40 degrees 
of freedom and forcing F — 8. Nevertheless, some newly obtained empirical closure 
equations expressing the average energy and the average momentum of the unperturbed 
system as simple power laws of F (with no dependence on N) have allowed to extend 
our results to the entire class of chaotic Lorenz 96 models. 

In this paper we have only used the KK relations in the most simplistic framework, 
i.e., computing the KK transforms and evaluating their agreement with the original 
data. Actually, several more sophisticated analysis techniques are available, such as 
recursive self-consistent algorithms, where the measured data are taken as first guess, 
exploiting the fact that multiple applications of KK relations, combined with sum rules, 
automatically filter our the noise and remove most of the spurious signal [40] . 

We believe that the proposed approach, which we may dub as spectroscopy of the 
climate system, may constitute a mathematically rigorous and practically very effective 
way to approach the problem of evaluating climate sensitivity and climate change from 
a radically new perspective. In this regard, we have proposed a rigorous way to compute 
the surface temperature response to changes in the the CO2 concentration at all time 
scales using only a specific set of simulations, and taking advantage of the theoretical 
results presented here. We underline that our approach takes into account all the 
(linear and nonlinear) feedbacks of the system, as they are included in the definition 
of the Green function. This, at a very practical level, is the great advantage of using 
Ruelle's formulas. 

At a more basic level, whereas considering more complex models requires heavier 
computational resources, the modest cost of the present set of simulations suggests 
that, at least for global or regional climatic observables, it is feasible to test the the- 
ory discussed here for simplified yet Earth-like climate models without resorting to 
top-notch computing facilities. Moreover, while in this paper we have computed the 
susceptibilities using, on purpose, a very cumbersome method, more efficient strategies 
can be devised, at least when the linear regime of the response is considered. Apart 
from the practical example given for the case of the impact of the CO2 concentration, 
these include studying the response of the system to S(t) like perturbations, which 
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gives directly the Green function of the system, and including in the forcing various 
monochromatic signals. Of course, in all cases, a Monte Carlo approach is needed in 
order to sample effectively the attractor of the unperturbed system in terms of the 
initial conditions of the simulations. 

These results pave the way for future investigations aimed at improving and extend- 
ing the theoretical framework presented here, at finding results of general applicability 
in the context of the modelling of geophysical fluid dynamics, and, finally at answering 
specific questions of relevance for climate dynamics. In this paper we have analyzed the 
simple case of the linear response, but, as discussed in [25, 17, 43], we have the algo- 
rithm to compute higher order terms, so that the treatment of the nonlinear response 
in entirely feasible. 

In the first direction, we foresee the possibility of writing out explicitly the linear 
susceptibility of a general observable by projecting the perturbation onto the unstable, 
neutral and stable manifolds and analyzing separately the contributions to the total 
response. This will probably require the adoption of adjoint techniques. Moreover, we 
will be testing the radius of convergence of the Ruelle response theory is some specific 
examples. 

Along the second direction, we propose to study the impact of stochastic forcing 
to deterministic chaotic models by treating the (additive or multiplicative) noise as 
a perturbation to be analyzed using the linear and nonlinear Ruelle response theory 
and related spectral methods. Moreover, we shall look into the spectral peaks of the 
susceptibilities and try to understand how the amplification of the response is related 
to resonances of the system and to the activation of positive feedbacks. 

Along the third direction, we envision the analysis of the impact of topography on 
the statistical properties of the circulation in a quasi-geostrophic setting, thus extending 
in a climatic perspective what presented in [64]. Moreover, we will tackle in an idealized 
setting the problem of computing the response of the storm track to changes in the 
surface temperature [65]. Moreover, we will try to compute along the way discussed 
here the climate response to changes in CO2 and solar irradiance using simplified but 
rather valuable climate models like PLASIM [66]. 

Finally, we would like to remark that the theory and the practical recipes proposed 
here could be of direct interest for all projects aimed at auditing climate models' per- 
formances and at studying practical problems connected to climate change, such as 
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PCMDI/CMIP3 (http://www-pcmdi.llnl.gov/), distributed computing initiatives 
such as climateprediction . net, and the new project PCMDI/CMIP5 (http : / /cmip- 
pcmdi . Ilnl.gov/cmip5/), which will provide a crucial input for the Fifth Assessment 
Report of the IPCC. 
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Figure 1: Linear susceptibility of intensive energy of the system E/N with respect to 
the global perturbation with X\ — 1 Vi. The real and the imaginary parts are depicted 
in Panels a) and b), respectively. The measured and extrapolated values are shown 
in red and black lines, respectively. The result of the Kramers-Kronig inversion done 
with the measured and with with the extrapolated data are shown in blue and magenta 
lines, respectively. 
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Figure 2: Linear susceptibility of intensive momentum of the system M/N with respect 
to the global perturbation with Xi = 1 Vi. The real and the imaginary parts are 
depicted in Panels a) and b), respectively. The measured and extrapolated values are 
shown in red and black lines, respectively. The result of the Kramers-Kronig inversion 
done with the measured and with with the extrapolated data are shown in blue and 
magenta lines, respectively. 
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Figure 3: Sum Rules for the real part of the linear susceptibility of intensive energy 
E/N (black line) and momentum M/N (blue line) of the system with respect to the 
global perturbation with X{ = 1 Vz. The theoretical values are indicated in the figure. 
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Figure 4: Green functions describing the time-dependent response of the observable 
of intensive energy E/N (black line) and momentum M/N (blue line) with respect to 
the global perturbation with Xi = 1 Vi. For both observables, the short time behavior 
of the Green function estimated from the asymptotic behavior of the corresponding 
susceptibility is shown in figure. 
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Figure 5: Comparison between the linear susceptibility of the intensive energy for the 
global perturbation and of the total energy for the local perturbation. The signals are 
the same except for the different level of noise. 
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Figure 6: Linear susceptibility of the energy at the grid point where the local pertur- 
bation is applied. The real and the imaginary parts are depicted in Panels a) and b), 
respectively. The measured and extrapolated values are shown in red and black lines, 
respectively. The result of the Kramers-Kronig inversion done with the measured and 
with with the extrapolated data are shown in blue and magenta lines, respectively. 
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Figure 7: Linear susceptibility of the momentum at the grid point where the local 
perturbation is applied. The real and the imaginary parts are depicted in Panels a) 
and b), respectively. The measured and extrapolated values are shown in red and black 
lines, respectively. The result of the Kramers-Kronig inversion done with the measured 
and with with the extrapolated data are shown in blue and magenta lines, respectively. 
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Figure 8: Sum rules of the real part of the linear susceptibilities indicated in the legend. 
The theoretical values are indicated in the figure. 
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Figure 9: Properties of the square of the linear susceptibility x £ n- The rea ^ anc ^ 
imaginary parts of [x^jv] 2 with their KK transforms are depicted in the first panel, the 
vanishing sum rule for the real part and the non-vanishing sum rule for the imaginary 
part are depicted in the second panel. 
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